Reconstructing image data from photon-counting measurements may for example include computing or estimating a vector a=(ajn) of basis coefficients of basis j=1, . . . , Me and image pixel n=1, . . . , Mp. based on a vector of photon counts m=(mik) in energy bin i=1, . . . , Me and detector element k=1, . . . , Md.
The photon counts may for example be obtained by direct readout pf a photon-counting x-ray detector. In another example, the photon counts may be the output of a post processing scheme operating on the counts obtained by readout from a photon-counting detector. The post-processing scheme may for example involve summation, filtering, averaging, and application of correction factors or correction terms.
In another example, the image data may consist of a single image, consisting of for example one of the basis coefficients or a combination thereof. In yet another example, the image data may consists of coefficients in a non-pixelized representation of an image. For example, the elements of a may be Fourier coefficients, or wavelet coefficients, or coefficients in a representation of the image as a sum of blobs.
The vector a is, in a typical case, obtained through the optimization of a function Φ(a, m) of a and m. A function which is optimized to find a will be referred to as a functional.
Such a functional, which is a function of measured data in the projection domain, i.e. image data relating to the transmission through or photon count values after an object, will be referred to as a projection-based functional. A projection-based functional may be used for image-based material decomposition or reconstruction. In another example, a projection-based functional may be a function of line integrals of the linear attenuation coefficient or of basis coefficients along a projection line.
A functional which is instead a function of measured data or data reconstructed from measured data, in the image domain, is called an image-based functional. Such a functional can be used to perform image-based material decomposition or reconstruction.
A functional Φ(a, m) may incorporate prior information about the imaged object. For example, this prior information may be provided in the form of an edge-preserving regularizer, which penalizes rapid variations in the image. In another example, the prior information may be provided in the form of a discrepancy term penalizing a difference between the reconstructed image and a prior image.
As an example, in the absence of pulse pileup, i.e. at low photon fluence rates, the number of registered photons in energy bin k and detector element k may be modelled as Poisson distributed with mean
                              λ          ik                =                                            ∫              0              kVp                        ⁢                                                            N                  0                                ⁡                                  (                  E                  )                                            ⁢                                                S                  i                                ⁡                                  (                  E                  )                                            ⁢                              exp                ⁡                                  (                                      -                                                                  ∑                                                  j                          =                          1                                                                          M                          b                                                                    ⁢                                                                                          ⁢                                                                                                    A                            jk                                                    ⁡                                                      (                                                          x                              ,                              y                                                        )                                                                          ⁢                                                                              μ                            j                                                    ⁡                                                      (                            E                            )                                                                                                                                )                                            ⁢              dE                                +                      s            ik                                              (        1        )            with Ajk=∫lk ajdl integrated along projection ray k. Here, N0(E) is the incident spectrum, Si(E) is a weight function modeling the sensitivity of the energy bin to different incident energy levels and sik is the number of expected scatter counts from the imaged object.
Eq. (1) does not model cross-talk between detector pixels, e.g. due to charge sharing, K-fluorescence or internal Compton scatter, but this can be included in the model by applying a linear operator B to the vector λ of expected counts: λc=B λ. For example, B may be represented by a sparse matrix with elements that are nonzero only for nearest-neighboring pixels. In another example, the elements of B may decrease with increasing distance between pixels. In yet another example, B may be obtained from a Monte Carlo simulation of a pencil beam of photons impinging on one detector element and scattered into neighboring detector elements.
Furthermore, Eq. (1) does not include pile-up. The effect of pile-up can for example be modelled by letting the registered counts in one projection line be a function of the incident photons in all different energy bins in the same projection line: λikp=f (λ1kc, . . . , λMe,kc) where λikp is the number of registered count after pile-up in energy bin i and detector element k. For example, λikp may be given by a paralyzable model or a non-paralyzable detector model.
Modern CT reconstruction algorithms typically generate the reconstructed image as a maximum a posteriori (MAP) estimate of the image given the measured data. The MAP estimate may build on a complete model of the relationship between the registered counts and the image values, or it may build on a simplified model of the relationship in order to simplify the optimization algorithm. Compared to the analytical reconstruction algorithms used previously, MAP reconstruction reduces noise and allows correction for detrimental effects such as scatter and optical blur. For energy-resolving, photon-counting CT, the MAP estimator of the vector of basis coefficients in the image a is, in a typical case:
                              a          ^                =                                            argmin              a                        ⁢                                          ∑                                  i                  =                  1                                                  M                  e                                            ⁢                                                          ⁢                                                ∑                                      k                    =                    1                                                        M                    d                                                  ⁢                                                                  ⁢                                  (                                                                                    λ                        ik                        p                                            ⁡                                              (                        a                        )                                                              -                                                                  m                        ik                                            ⁢                                                                        λ                          ik                          p                                                ⁡                                                  (                          a                          )                                                                                                      )                                                              +                      R            ⁡                          (              a              )                                                          (        2        )            where λikp and mik are the number of expected and registered counts, respectively, in energy bin i=1, . . . , Me and detector element k=1, . . . , Md. λikp includes the effect of cross-talk and pileup is obtained from (1) with the blur operator B and the pileup function ƒ applied. R(a) is an edge-preserving regularizer, which penalizes differences between neighboring detector elements. The expression that is optimized in (2) will be referred to as an MAP functional.
If the regularization term is not included, (2) becomes a maximum likelihood (ML) estimator, and the functional to be minimized is called a maximum likelihood functional
Since (2) is difficult to solve fast enough, it is common practice to replace it with a simplified penalized weighted least squares estimator. This estimator is given by minimization of a penalized weighted least squares functional, for example:
                              a          ^                =                                            argmin              a                        ⁢                                          ∑                                  j                  =                  1                                                  M                  b                                            ⁢                                                          ⁢                                                ∑                                      k                    =                    1                                                        M                    d                                                  ⁢                                                                  ⁢                                                                            (                                              Ta                        -                                                                              A                            ^                                                    jk                                                                    )                                        2                                                                              σ                      2                                        ⁡                                          (                                                                        A                          ^                                                jk                                            )                                                                                                    +                      R            ⁡                          (              a              )                                                          (        3        )            
Here, T denotes the forward ray transform operator and Âjk is an estimate of the line integral Ajk=∫lk ajdl along projection ray k. σ2(Âjk) is the variance of Âjk. Âjk can be obtained from the measured counts mik for each individual detector element using maximum likelihood estimation or a look-up table. Eq. (3) can be computed quickly using e.g. the iterative coordinate descent (ICD) method or the separable quadratic surrogates (SQS) method, but it gives inferior image quality compared to (2) since it builds on a simplified noise model and ignores detector cross-talk and object scatter. In particular, (3) is based on modelling the noise as Gaussian instead of Poisson, which is an approximation. There is thus a need for an algorithm that combines good image quality with fast computation time.
The publication “Multi-Material Decomposition Using Statistical Image Reconstruction for Spectral CT” by Y. Long and J. Fessler, IEEE Trans. Med. Imag. 33, pp. 1614-1626 (2014), relates to an iterative reconstruction method for penalized-likelihood reconstruction of material basis images. This method is initialized by a set of basis images obtained through an image-domain material decomposition method.
U.S. Pat. No. 5,390,258 relates to a method of acquiring an image from an object, wherein a set of training images is used to generate a convergent series expansion, and wherein the measured signals are used to generate a truncated series expansion of an image of the object.
U.S. Pat. No. 6,754,298 relates to a method based on a statistical model for reconstructing images from transmission measurements with high energy diversity.
U.S. Pat. No. 7,551,708 relates to a method of reconstructing material decomposed images from data from energy discriminating computed tomography detectors using the iterative coordinate descent (ICD) algorithm.
U.S. Pat. No. 9,165,384 relates to a method of image reconstruction which reconstructs a plurality of final component images of an object based on spectral projection data, wherein intermediate images are used in the reconstruction algorithm and wherein correlations between these intermediate images are taken into account in the algorithm.
U.S. Pat. No. 8,923,583 relates to a tomographic reconstruction method wherein material component images are reconstructed by optimizing a joint likelihood functional, which includes information on correlations between the component sinograms.
US patent application 2016120493A1 relates to an x-ray CT image processing method wherein a joint posterior distribution, based on a prior probability distribution, is used to estimate an x-ray absorption coefficient from measurements with different wavelengths.
U.S. Pat. No. 8,929,508 relates to a method of computing line integrals of basis coefficients through an object from x-ray photon transmission measurements by computing a first approximation to the line integrals and combining the first approximation with a correction computed from a calibration phantom.
U.S. Pat. No. 6,907,102 relates to a method of image iterative reconstruction wherein a cross-section reconstruction vector approximately matching the projection data is determined using a computed tomography model.
U.S. Pat. No. 7,885,371 relates to a method of tomographic image reconstruction wherein a first reconstruction method, which converges faster on low spatial frequencies than on high spatial frequencies, is followed by a second reconstruction method which converges faster on high spatial frequencies than on low spatial frequencies.
U.S. Pat. No. 9,508,163 relates to a method of iterative tomographic reconstruction wherein each iteration of an outer loop includes iterative processing of an inner loop.
U.S. Pat. No. 9,020,230 relates to a reconstruction method employing an outer iteration loop and an inner iteration loop, wherein the inner iteration loop calculates a preconditioner used by the outer loop.
U.S. Pat. No. 6,256,367 relates to a method of correcting for artifacts due to scatter CT images by Monte Carlo simulating photon scatter and subtracting the simulated photon energy from the measured projection data.
There is, however, an ongoing need for a fast reconstruction algorithm that gives good image quality.